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ABSTRACT 


^-This  thesis  analyses  data  o-f  and  builds  a  simulation 
model  for  the  track  o-f  an  underwater  vehicle  as  perceived 
by  a  test  range  o-f  three  dimensional  short  baseline  sonar 
arrays.  In  this  way  many  random  replications  of  track 
become  available  quickly  and  inexpensively.  These 
simulations  support  a  larger  project  whose  object  is  to 
monitor  the  performance  of  the  test  range  and  provide  clues 
for  troubleshooting  problems.  In  particular,  joint  values 
of  sensor  array  estimated  displacement  and  reorientation 
corrections  are  generated  and  their  error  distribution  is 


quantified 
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I.  INTRODUCTION 


The  Naval  Undersea  Weapons  Engineering  Station  (NUWES) 
operates  nine  test  ranges  that  serve  a  variety  of  purposes 
for  testing  and  validating  the  Navy's  current  and  future 
weapons  systems.  This  thesis  is  in  support  of  a  larger 
project  whose  goal  is  to  monitor  the  performance  of  the 
short  baseline  ranges  (ie.  those  ranges  in  which  each  array 
produces  a  three  dimensional  track  of  moving  vehicles). 
More  specifically,  it  deals  with  the  development  of  a 
computer  simulation  which  generates  realistic  random 
replications  of  an  underwater  vehicle's  track  (ie.  vehicle 
position  perceived  by  a  sensor  array  on  the  range  at  a 
number  of  equally  spaced  time  points).  Such  simulations 
serve  to  provide  information  about  the  inherent  variability 
in  the  output  of  the  tracking  range,  especially  for 
assessing  the  calibration  error.  Its  use  can  provide  an 
improved  understanding  of  the  processes  involved  and  a 
reduction  in  the  frequency  of  range  shutdown  for  the  purpose 
of  resurveying  the  remote  sensor  locations. 

When  the  simulated  replicate  tracks  are  passed  to  the 
program  kEYMAIN  (a  project  FORTRAN  program  that  estimates 
displacement  and  orientation  corrections  to  the  range  remote 
sensor  arrays) ,  the  output  is  used  to  generate  bivariate 


scatter  plots  of  these  corrections.  Although  extensive  work 


of  this  type  was  not  possibla  in  ths  currant  thesis,  the 
limited  results  indicate  a  surprising  amount  of  variability 
and  suggest  that  the  nature  and  extent  of  variability  can 
change  noticeably  with  relatively  minor  changes  in  the 

localized  conditions.  Considerable  testing  using  this 
simulation  tool  is  clearly  indicated. 

The  research  reported  here  involves  three  general 
activities: 

(a)  Analysis  of  real  underwater  track  data  to  learn  their 
important  behavioral  characteri sti cs. 

(b)  Simulation  model  formulation  and  construction. 

(c)  Development  and  programming  of  algorithms  to  merge 

with  the  existing  project  programs  to  produce  any 

number  of  random  replications  with  correction 
estimates  for  each  simulated  track. 

The  organization  of  the  thesis  is  as  follows: 

Section  II  contains  explicit  background  material  and 
provides  a  framework  for  the  research  and  shows  how  it 
relates  to  the  larger  project.  Section  III  explains  the 
data  analysis  procedures  and  results.  The  simulation  model 
is  developed  in  Section  IV  and  the  interfacing  with  the 

overall  project  programs  is  described  in  Section  V.  Results 
and  conclusions  are  detailed  in  Section  VI,  with  areas  for 
future  work  listed  in  Section  VII. 

A  number  of  appendices  are  included  to  provide  the 

detailed  support  too  extensive  to  be  included  in  the  main 
body  of  the  paper.  They  are  referenced  in  the  appropriate 


sections. 


Additional ly. 


to  speed  the  completion  and 


availability  of  the  KEYMAIN  program,  the  author  wrote  two 
subroutines  (CONECT  and  REDUCE)  to  be  included  in  that 
package.  Appendices  E  through  G  contain  the  development  of 
these  subroutines  and  the  FORTRAN  code  listings. 


II.  BACKGROUND 


Consider  a  three  dimensional  underwater  tracking  array 
composed  of  four  hydrophones  arranged  to  define  a  local 
Cartesian  coordinate  system  (see  Figure  1).  Such  an  array 
is  capable  of  tracking  a  target  downrange,  crossrange  and  in 
depth  in  its  local  coordinate  system;  it  is  termed  a  short 
baseline  array  because  of  the  short  distance  of  the  X.  Y, 
and  Z  hydrophones  from  the  '‘corner"  hydrophone  (30  feet) 
that  is  used  to  gather  the  three  dimensional  tracking 
i  nf  or  mat  i  on . 

The  arrays  take  "fixes"  of  target  position.  A  fix  is 
defined  by  a  bearing  (azimuth  and  elevation)  and  distance 
from  the  array  to  the  target.  Tracked  targets  on  the  range 
are  fitted  with  an  acoustical  device,  called  a  "pinger " , 
that  emits  pulses  of  sound,  or  "pings",  at  a  specific 
frequency  at  precisely  timed,  regular  intervals.  The 
elapsed  time  from  the  generation  of  the  ping  at  the  target 
to  the  reception  of  the  ping  at  the  array  establishes  a 
distance  from  the  target  to  the  array.  Because  of  the 
distance  that  separates  the  individual  hydrophones  of  the 
array,  the  ping  arrives  at  each  hydrophone  at  a  slightly 
different  time.  This  time  difference  can  be  resolved  into  a 
bearing  from  the  array  to  target.  (Actually,  these  fixes 
are  "apparent"  fixes.  They  are  used  to  initialize  a  sound 

1  1 


ray  tracing  algorithm  that  leads  to  the  "real"  fix.)  A 
series  of  such  fixes  establishes  a  target  track  in  that 
particular  array's  own  local  Cartesian  coordinate  system. 
If  the  precise  position  ni  the  array  is  known,  its  local 
track  information  can  be  translated  to  one  common  range 
coordinate  system. 

Position  of  an  individual  array  in  the  range  Cartesian 
coordinate  system  is  described  not  only  in  terms  of  its 
downrange,  crossrange  and  depth  (X,Y  and  Z)  coordinates 
(termed  location),  but  also  with  respect  to  X-tilt,  Y-tilt 
and  Z-rotation  (commonly  called  roll,  pitch  and  yaw)  angles 
from  the  range  coordinate  system  axes  (termed  orientation). 
Both  sets  of  measures  are  needed  to  translate  accurately 
from  local  to  range  coordinates. 

Typically,  a  range  is  composed  of  many  individual 
arrays.  Nanoose  Range,  for  example,  has  24  arrays,  while 
Dabob  Bay  has  7.  The  arrays  are  arranged  in  such  a  way  that 
the  array  coverages  overlap  one  another  to  provide 
continuous  tracking  on  a  target  vehicle.  Such  overlapping 
areas  are  called  crossover  regions  (see  Figure  2),  and 
produce  two  sets  of  track  on  the  same  target  for  the  same 
time  period. 

Ideally,  the  correspondi ng  points,  or  fixes  described 
earlier,  from  each  array  in  a  crossover  region  should 
translate  to  identical  points  on  the  range  coordinate 
system.  In  practice,  this  seldom  happens. 


SENSOR  ARRAY  CROSSOVER  REGIONS 


. NOT  TRACKED 

—  SINGLE  COVERAGE 

—  DOUBLE  COVERAGE 


Figure  2  :  Crossover  Regions 


Three  major  sources  of  this  variability  in  crossover 
track  data  ares 

(a)  Slippage  of  the  sensor  arrays  from  their  assumed 
positions  in  the  range  coordinate  system. 

(b)  Time  synchronization  and  instrumentation  problems. 

(c)  Inhomogeneities  and  temporal  variability  in  the 
water  column. 

Although  the  reasons  for  slippage  of  the  arrays  are 
speculative,  the  fact  that  slippage  occurs  is  evidenced  by 
the  change  in  sensor  locations  after  a  range  resurvey  is 
performed.  The  question  of  discriminating  timing  errors, 
item  (b) ,  from  slippage  errors  is  on  a  future  agenda. 
Investigation  done  by  Main  CRef.  11  provides  us  with 
methodology  for  treating  the  random  components  of  these 
errors.  Item  <c>  falls  into  a  broader  category  which  may 
require  extensive  investigation.  It  may  also  provoke  a 
reassessment  of  the  assumed  uniform  horizontal  quality  of 
the  range  water  column.  It  seems  wise  to  account  for 
slippage  and  instrumentation  problems  first.  Our  work  is 
confined  to  (a) . 

The  present  research  is  in  direct  support  of  Professor 
Robert  R.  Read  of  the  Naval  Postgraduate  School  who  has  been 
working  on  the  slippage  question.  He  has  produced  a  FORTRAN 
program  (KEVMAIN)  which  takes  the  crossover  data  and 
estimates  array  position  corrections  using  a  non  linear 
least  squares  algorithm.  The  surfaces  that  enter  into  this 
least  squares  optimization  function  are  rather  flat  and  the 
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values  do  not  change  much  a«  the  location  corractions  are 
variad.  Replicated  data  is  needed  to  quantify  the  extent  to 
which  the  estimated  corrections  are  within  tha  range  o-f 


natural  variability.  The  simulation  model  in  this  thesis 
provides  the  tool  by  which  this  variability  can  be 
quantified.  With  the  simulation,  the  needed  replicated  data 
can  be  quickly  and  easily  generated,  and  scatterplots  of 
displacement  (change  in  location)  versus  rotation  (change  in 
orientation)  can  be  made  to  investigate  the  inherent 
variability  in  the  correction  estimates. 


A.  GENERATING  RESIDUALS  FROM  REAL  TRACK  DATA 

The  real  track  data  used  for  this  study  was  taken  from 

Nanoose  range  in  September ,  19B2  (Figure  3).  It  was 

« 

supplied  as  an  N  x  7  matrix,  N  representing  the  number  of 
data  points  in  the  crossover  data  set  and  the  columns 
representing  the  following: 

Column  1  Point  count.  Every  pulse  emitted  by  the 

pinger  is  assigned  a  sequential  number 
beginning  at  the  start  of  the  tracking  run. 
Therefore,  this  column's  value  increases  by 
one  each  row  unless  a  data  point  is  missing, 
which  would  be  indicated  by  a  sequential 
omission  in  this  first  column. 

Columns  2-4  The  X,Y,  and  Z  coordinates  of  the  target 
for  the  column  one  point  count  as  determined 
by  the  first  sensor  in  the  crossover  pair. 

Columns  5-7  The  X,Y,  and  Z  coordinates  of  the  target  for 
the  column  one  point  count  as  determined  by 
the  second  sensor  in  the  crossover  pair. 

There  were  six  such  data  sets  used.  These  particular 
sets  were  chosen  because  each  formed  a  straight  line  in  the 
range  space  and  so  the  best  straight  line  path  of  the  target 
could  be  estimated  from  the  data.  (In  contrast,  estimating 
the  best  curved  path  from  curved  track  data  would  have  been 
far  more  challenging,  but  with  little  or  no  gain  in  the 
examination  of  residuals.) 

The  idea  was  to  take  the  track  data,  fit  the  best 
straight  line  possible  through  the  data,  then  examine  the 


residuals  formed  by  subtracting  th®  fitt®d  straight  lin® 
point  sequences  from  the  original  data. 

The  commented  FORTRAN  program  written  to  generate  the 
residuals  is  included  as  Appendix  A,  and  its  mechanics  are 
discussed  below. 

Since  the  data  was  given  in  a  single  matrix  and 
residuals  were  desired  for  each  sensor,  the  data  was  first 
separated  into  two  tracks,  each  an  N  x  3  matrix;  call  them 
TR1  and  TR2.  Next,  a  3  x  3  covariance  matrix  for  each  track 
was  computed  by  the  formula 

N  __ 

COV  *  C  5H  (TR(i)-TR)'  (TR(i)-TR)3  /  N-l 

i  =  1 

where 

fR(i)  =3  component  row  vector  of  the  X,Y,Z  coordinates  of 
the  i data  point 

TR  =  3  component  row  vector  of  the  average  of  the  N 
( X , Y , Z )  coordinates  for  that  track  data 
N  =  number  of  data  points 

The  best  straight  line  estimate  of  the  target  path  is 
the  straiqht  line  such  that  the  sum  of  the  distances  of  each 
data  point  to  the  line  is  a  minimum.  This  implies  that  the 
distance  from  each  point  to  the  line  is  a  minimum  -  that  is, 
each  point  should  be  projected  onto  the  line  orthogonally. 
The  method  employed  to  accomplish  this  orthogonal  regression 
was  that  of  principle  components.  Principle  components 
requires  that  the  eigenvector  associated  with  the  largest 
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eigenvalue  of  the  covarianc*  matrix  be  idanti-fiad.  It  was 
computationally  convenient  to  make  a  3  x  3  matrix  o-f  the 


eigenvectors,  with  the  -first  column  being  the  eigenvector 
associated  with  the  largest  eigenvalue,  the  second  column, 
the  eigenvector  associated  with  the  second  largest 
eigenvalue  and  third  column,  the  eigenvector  associated  with 
the  smallest  eigenvalue.  This  done,  the  -following  matrix 
multiplication 

PROJECTION  -  <TR  -  AVE)  x  EIGENVECTORS 

where 


TR  -  N  x  3  track  data  matrix 
AVE  “  N  x  3  matrix  made  o-f  N  identical  rows  o-f 
the  X,Y,  and  Z  averages  o-f  the  track  data 
EIGENVECTORS  *  3  X  3  eigenvector  matrix  described  above 
yields  the  N  x  3  matrix  PROJECTION,  whose  entries  are  the 
orthogonal  projections  o-f  the  track  data  points  onto  the 
axes  of  a  new  3  dimensional  Cartesian  coordinate  system 
whose  origin  is  located  at  TR  (the  <X,Y,Z>  averages  for  the 
data  set)  and  whose  axes  are  rotated  such  that  the  new  X 
axis  is  the  best  straight  line  that  describes  target  path. 
For  example,  the  first  row  of  PROJECTION  is  a  three 
component  row  vector;  call  the  components  <pl,p2,p3).  Then 
the  point  (pi, 0,0)  is  the  projection  of  the  first  track  data 
point  onto  the  new  X  axis  (known  to  be  the  best  straight 
line  that  describes  the  target  path),  the  point  (0,p2,0)  is 
the  projection  of  the  first  track  data  point  onto  the  new  Y 
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axis,  and  tha  point  <0,0, p3)  is  the  projection  o-f  the  -first 
track  data  point  onto  the  new  Z  axis.  Since  we  are  only 
interested  in  the  projections  onto  the  new  X  axis,  we  can 
replace  the  second  and  third  columns  o-f  PROJECTION  with 
zeros  and  have  the  coordinates  o-f  the  best  straight  line 
path  o-f  the  target.  These  coordinates  are  in  the  PROJECTION 
coordinate  system,  and  must  be  translated  back  into  the 
range  coordinate  system.  Using  the  matrix  PROJECTION  (with 
last  two  columns  ■  0)  and  performing  the  following 

mul tipi i cat ion 

OLD  -  EIGENVECTORS  x  PROJECTION' 
yields  the  3  x  N  matrix  of  the  track  data  points  in  the 
range  coordinate  system.  Note  that  OLD  must  be  transposed 
to  get  the  N  x  3  format  desired.  Since  the  mean  had  been 
subtracted  off  before  computing  the  PROJECTION  matrix,  to 
get  the  data  paints  back  to  the  proper  position  requires 
that  the  means  be  added  back  in.  The  matrix  addition 

OLD '  +  AVE 

yields  the  N  x  3  matrix  of  the  orthogonal  projections  of  the 
data  onto  the  straight  line  path  of  the  target  in  the 
original  range  coordinate  system. 

Actually,  after  having  computed  the  PROJECTION  matrix, 
there  is  yet  another  step  to  take  before  translating  back 
into  the  range  coordinate  system.  Based  on  an  assumption  of 
constant  speed  for  the  target,  the  track  data  paints  should 
be  equally  spaced  since  the  pulses  are  emitted  at  regular. 


precisely  timed  intervals.  This  could  be  accomplished  by 
using  an  interval  equal  to  the  total  distance  divided  by  the 
number  of  data  points,  but  if  data  points  are  missing  in  the 
track  data  set  (and  each  set  of  track  data  was  missing 
several  points),  then  this  method  introduces  an  error.  The 
method  chosen  should  show  that,  if  a  single  data  point  is 
missing,  the  gap  between  the  consecutive  points  on  the  data 
track  should  be  twice  as  long  as  if  none  were  missing.  The 
first  column  (  the  point  counts)  of  the  N  x  7  track  data 
matrix  contains  the  information  on  missing  points. 
Performing  a  simple  least  squares  regression  of  the  point 
counts  onto  the  first  column  of  the  PROJECTION 
matrix transl ates  this  information  into  the  best  straight 
line  point  sequence  for  target  motion.  Now  when  the 
projections  onto  the  first  principle  component  are 
translated  back  onto  the  range  coordinate  system,  the 
result  is  truly  the  best  set  of  track  points  attainable 
using  all  the  information  contained  in  the  track  data  set. 

To  obtain  the  residuals,  the  straight  line  estimated 
track  points  are  subtracted,  point  for  point,  from  the 
original  track  data.  It  was  important  to  discover  the 
distribution  of  the  residuals  for  each  track  segment  because 
to  simulate  the  track  segments  later,  residuals  were 
simulated  and  then  added  to  the  straight  line.  Knowing  the 
character  of  the  residuals  allowed  the  generation  of  very 


realistic  track  data  for  the  simulation  model. 


B.  ANALYZING  TRACK  DATA  RESIDUALS 


Each  of  the  six  track  segments  produced  two  sets  o-f 
residuals  -  one  set  for  each  of  the  two  sensors  of  the 
crossover  pair.  Each  set  was  further  divided  into  X,  Y,  and 
Z  components.  There  were,  therefore,  36  sets  of  residuals  to 
be  examined  for  distributional  characteri sties.  Graphical 
analysis  was  performed  on  each  set  of  residuals  (histograms, 
cumulative  density  plots  and  QQ  plots)  and,  when,  from  that 
analysis,  a  distribution  for  the  residuals  could  be 
determined,  formal  statistical  tests  were  performed  to 
verify  that  the  best  distribution  was  chosen.  The  analysis 
showed  the  residuals  to  be  normally  distributed.  The  best 
and  worst  case  graphical  and  analytical  results  are  included 
in  Appendix  B.  Note  that  there  were  some  very  good  fits  to 
a  normal  density  <pp.  66-71)  with  statistical  significances 
for  the  Chi  Squared  and  Kolmogorov-  Smirnov  goodness  of  fit 
tests  well  above  a  very  conservative  .35  in  all  but  one 
case.  Even  in  the  data  that  most  poorly  resembled  a  normal 
density  <pp.  72-77) ,  the  Kolmogorov  -  Smirnov  goodness  of 
fit  significance  level  never  fell  below  .25.  From  this 
analysis  it  was  concluded  that  a  normal  density  for  the 
residuals  was  accurate  and  appropriate  for  the  simulation 
model . 

C.  TIME  SERIES  ANALYSIS 

From  a  simulation  point  of  view,  it  is  very  convenient 
to  assume  independence  of  the  residuals  between  successive 
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points  in  the  original  track.  If  the  assumption  is  not 
made,  then  one  must  resort  either  to  using  a  point  count 
interval  of  some  integer  value  greater  than  one  for  making 
the  simulated  track  from  the  original,  or  build  an 
autoregressive  model. 

Concern  for  this  correlation  of  the  data  over  time  led 
to  a  time  series  analysis  o-f  the  residuals.  This  was  done 
in  two  ways  for  each  track  from  a  single  sensor.  Recall 
that  for  each  data  point,  residuals  were  generated  in  the  X, 
V,  and  Z  directions  of  the  range  coordinate  system.  Each  of 
these  components  was  subjected  to  a  time  series  analysis  to 
determine  whether  there  was  a  dependence  in  the  errors  in 
any  single  direction  over  time.  Then,  the  distance  of  the 
data  point  from  the  straight  line  target  path,  given  by 

SQRTC  (RES  2)  +  (RES  2)  +  (RES  2>  1 
x  y  z 

was  examined  to  determine  if  the  magnitude  of  the  error  of 
one  point  was  correlated  with  the  magnitude  of  the  error  of 
the  next  point.  The  time  series  analysis  examines  the 
dependence  from  point  to  point,  on  every  other  point,  on 
every  third  point,  and  so  on,  so  that  the  interval  at  which 
one  can  assume  independence  (indicated  by  an  autocovariance 
value  of  zero)  can  be  determined.  Appendix  C  contains  the 
autocorrel at i on  graphs  for  the  6  data  sets;  first  the  error 
magnitudes  are  examined  for  both  sensors  of  a  data  set, 
followed  by  the  individual  error  analysis  in  the  X,  Y,  and  Z 


d i rec t i ons . 


Thera  was  no  clear  cut  interval  -for  which  independence 
seemed  to  emerge  in  the  data  sets.  There  is  instead  a 
random  pattern  of  insignificant  dependence  from  the 
beginning  of  the  autocorrel ati on  analysis,  leading  to  the 
conclusion  that  point  to  point  independence  was  appropriate 
for  the  simulation  model.  Attention  is  directed  to  the 
noise  in  the  autocorrel ation  graphs  with  no  distinct 
patterns  emerging,  and  correlation  magnitudes  smaller  than 
0.2  in  most  cases.  It  was  decided,  therefore,  that  any 
important  time  correlation  was  not  so  large  as  to  cause 
excessive  or  even  noticeable  error  in  the  simulated  tracks. 

D.  STUDY  OF  RESIDUALS 

The  study  of  the  residuals  yielded  important  interesting 
information,  summarized  in  Table  1.  The  first  column  of 
Table  1  is  the  standard  deviation  of  the  residuals  from  the 
left  sensor  in  the  downrange  (X) ,  crossrange  (Y)  and  depth 
(Z)  directions,  three  values  for  each  data  set.  Column  B 
gives  the  correspondi ng  information  for  the  right  sensor  of 
the  crossover  data  pair.  Columns  2,  3  and  4  for  each  data 
set  are  the  columns  of  a  3  x  3  correlation  matrix  for  the 
residuals  of  the  left  sensor,  columns  6  and  7, 
correspondi ng  information  for  the  right  sensor. 

First,  note  the  difference  between  the  left  and  right 
sensor  standard  deviations  of  residuals  in  all  three 
directions.  Although  there  are  some  very  good  comparisons 
(data  set  D2A1,  crossrange,  and  data  set  D5A ,  downrange) 
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TABLE  1 

STANDARD  DEVIATIONS  OF  RESIDUALS 
AND 

CORRELATION  COEFFICIENTS  OF  RESIDUALS 


SD  CORRELATION  CORRELATION  SD 

LEFT  LEFT  RIGHT  RIGHT 


DATA  SET  D2A1  ,  N  =  67 


.842 

1.000 

0.916 

-0.899 

1.000  -0.863 

0.787 

0. 

721 

.  141 

0.916 

1.000 

-0.728 

-0.863  1.000 

-0.455 

n 

m 

154 

.394 

-0.899 

-0.728 

1.000 

0.787  -0.455 

1 . 000 

r> 

063 

DATA  SET 

D2A2  ,  N  =  70 

.  484 

1 . 000 

0.859 

-0.787 

1.000  -0.360 

0 .040 

i . 

163 

.  396 

0.859 

1.000 

-0.759 

-0.360  1.000 

-0.547 

147 

.694 

-0.787 

-0.759 

1.000 

0.040  -0.547 

1 . 000 

i  . 

976 

DATA  SET 

D2B  ,  N  =  53 

.693 

1 . 000 

0.788 

-0.888 

1.000  -0.949 

0.902 

0. 

937 

.  OOO 

0.  783 

1.000 

-0.434 

-0.949  1.000 

-0. 723 

• 

350 

.  353 

-0.888 

-0. 434 

1 . 000 

0.902  -0.728 

1 . 000 

042 

DATA  SET 

D4  ,  N  =  93 

.  497 

1 . 000 

0.741 

-0. 904 

1.000  -0.681 

0 . 330 

0. 

502 

.955 

0.741 

1 . 000 

-0. 687 

-0.681  1.000 

-0. 146 

r~y 

4  b  7 

.537 

-0 . 904 

-0.687 

1 . 000 

<1 

■ 

o 

1 

o 

CD 

n 

■ 

o 

1 . 000 

i . 

843 

DATA  SET 

D5A  ,  M  =  82 

.  7S1 

1 . 000 

-0.930 

-0. 429 

1.000  -0.953 

0.168 

t 

i.  • 

T  O  ^ 

/  iju 

.  751 

-0.930 

1 . 000 

0.  152 

-0.953  1.000 

-0. 148 

<2  m 

Z-  0  " 
i— f  .„j 

.  333 

-0. 429 

0.  152 

1 . 000 

0.168  -0.148 

1 . 000 

1  . 

378 

DATA  SET 

D5B  ,  N  =  87 

.  386 

1 . 000 

-0 . 408 

-0. 450 

1.000  0.445 

0.  395 

0. 

71' 

.931 

-0. 403 

1 . 000 

-0.454 

0.445  1.000 

-0. 279 

1 . 

a  1  3 

.  397 

-0 . 450 

-0.454 

1 . 000 

0.395  -0.279 

1 . 000 

4. 

254 

2b 
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there  arc  also  soma  wide  disparitias  (data  set  D2B, 
crossrange  and  data  set  D5B,  crossrange) .  Using  the 
variance  ratio  test  described  in  Larson  CRef.  2:  pp.  4491,  a 
statistical  test  was  performed  to  determine  whether  the 
variances  (the  recorded  standard  deviations,  squared)  of  the 
left  and  right  sensor  arrays  for  downrange,  crossrange  and 
depth  were  the  same  for  a  given  data  set.  The  results  are 
given  in  Table  2.  Using  0.05  as  the  basis  for  rejection  of 
the  null  hypothesis  (H^:  the  two  variances  are  the  same),  8 
of  the  18,  or  44%,  of  the  individual  tests  fail.  With  a 
confidence  level  of  .95,  one  would  expect  roughly  1  failure 
in  the  18  tests.  Eight  failures  is  strong  evidence  that  the 
standard  deviation  figures  do  not  match  very  well. 

Recalling  the  test  data  from  Figure  3  (p.  18),  it  is 
seen  that  there  are  three  sensor  arrays  that  contributed  the 
data:  arrays  4,  5  and  6.  In  data  sets  D2A1 ,  D2A2  and  D2B, 
array  4  is  the  left  array  and  array  5  is  the  right  sensor 
array.  For  data  sets  D4,  D5A  and  D5B,  sensor  array  5  is  the 
left  array  and  array  6  is  the  right  array.  One  might  expect 
that  the  standard  deviations  in  any  single  direction  would 
be  the  same  for  a  given  sensor.  This  turned  out  not  to  be 
the  case.  Using  Bartlett's  test  for  the  equality  of  several 
variances  CRef.  3:  pp.  225-2273,  the  hypothesis  that  the 
variances  for  a  given  sensor  in  a  given  direction  are  equal 
was  tested.  The  results  are  given  in  Table  3.  The  first 


two  sections  of  the  table  compare  the  3  standard  deviation 


TABLE  2  TEST  OF  THE  HYPOTHESIS  THAT  THE 
STANDARD  DEVIATIONS  WITHIN  A  DATA  SET 


ARE 

EQUAL 

VARIANCE 

RATIO 

DEGREES 

OF  FREEDOM 

LEVEL  OF 
SIGNIFICANCE 

DATA  SET  D2A1 

DOWNRANGE 

1.363 

66 

.211 

CROSSRANGE 

.988 

66 

.935 

DEPTH 

1.363 

66 

.229 

DATA  SET  D2A2 

DOWNRANGE 

1.627 

69 

.045 

CROSSRANGE 

1.244 

69 

.364 

DEPTH 

1.858 

69 

.011 

DATA  SET  D2B 

DOWNRANGE 

.547 

52 

.031 

CROSSRANGE 

.357 

52 

.0003 

DEPTH 

.598 

52 

.067 

DATA  SET  D4 

DOWNRANGE 

.982 

92 

.884 

CROSSRANGE 

.628 

92 

.027 

DEPTH 

.692 

92 

.079 

DATA  SET  D5A 

* 

DOWNRANGE 

.995 

81 

.939 

CROSSRANGE 

.749 

81 

.  195 

DEPTH 


508 


81 


003 


TABLE  2 

(conti nu«d) 

DATA  SET  D5B 

DOWNRANGE 

1.527 

86 

.051 

CROSSRANGE 

3.300 

86 

0.0 

DEPTH 


638 


86 


038 


Table  3  TEST  OF  THE  HYPOTHESIS  THAT  THE 
STANDARD  DEVIATIONS  OF  RESIDUALS  OF 
A  PARTICULAR  SENSOR  ARE  EQUAL 


Chi  Squared  Random  Variables 
DEGREES  OF  FREEDOM  .99  QUANTILE  .999  QUANTILE 


2 

9.21 

13.82 

5 

15.09 

20.52 

DOWN 

RANGE 

CROSS 

RANGE 

DEGREES 

DEPTH  OF  FREEDOM 

SENSOR 

4  (LEFT) 

39.57 

2.02 

1.41 

2 

SENSOR 

5  (RIGHT) 

14.  75 

16.25 

13.96 

2 

SENSOR 

5  (LEFT) 

131.14 

101.33 

91.40 

2 

SENSOR 

6  (RIGHT) 

148.72 

175. IB 

84.31 

2 

SENSOR 

5  (LEFT  AND 
RIGHT) 

153.04 

150. B7 

105.51 

5 

LEFT  SENSORS 

171.45 

162.29 

94.48 

5 

RIGHT  SENSORS 

168.37 

236.97 

107. 15 

5 

p 


valuat  -for  each  imaor  in  a  singla  diraction.  Tha  tast 
statistic  is  distributad  as  a  Chi  Squarad  random  variabla 
with  2  dagraas  of  freedom.  Tha  .9?  and  .999  quantiles  for 
such  a  random  variabla  are  9.21  and  13.82.  In  only  two 
cases  of  tha  12  trials  would  tha  3  standard  deviations  be 
considered  statistically  tha  same. 

The  third  section  of  Table  3  is  included  because  sensor 
array  5  was  used  as  both  a  right  and  left  array.  Therefore, 
there  were  actually  six  values  for  downrange,  crossrange  and 
depth  for  that  particular  array.  Using  the  same  test  as 
before,  the  question  of  whether  the  six  values  were 
statistically  the  same  was  investigated.  The  tabulated 
statistic  for  the  3  cases  is  distributed  as  a  Chi  Squared 
random  variable  with  5  degrees  of  freedom.  The  .99  and  .999 
quantiles  for  such  a  ramdom  variable  are  15.09  and  20.52. 
In  every  case,  the  p-value  of  the  test  is  essentially  equal 
to  0.0.  The  results  are  highly  significant. 

Finally,  the  bottom  of  the  table  investigates  whether 
downrange,  crossrange  and  depth  residual  standard  deviations 
are  the  same  for  left  and  right  sensors  in  general. 
Employing  the  same  test  for  the  6  values  produced  the 
tabulated  statistics.  Since  the  test  statistic  is  again  Chi 
5quared  with  5  degrees  of  freedom,  the  level  of  significance 
in  every  case  is  essentially  0.0. 

A  third  question  of  interest  from  Table  1  is  whether 


the  correlation  coefficients  are  the  same  for  the  left  and 


right  sensors  in  a  data  sat.  This  was  invastigatad  in  tha 
following  mannar. 

First,  tha  normalizing,  inverse  hyperbolic  tangent 
transformation  CRef.  3:  p.  3651  was  applied  to  each  of  the 
off  diagonal  correlation  coef f iceients  in  the  correlation 
matrices.  (Note  that  the  matrices  are  symmetric,  so  there 
are  only  three  values  of  concern  for  each  matrix.)  This 
transf ormation  makes  each  of  the  values  normal  with  mean 


1  1  +  rho 

-  In  — — 

2  1  -  rho 

and  variance 


1 

N  -  3 

where  N  is  the  number  of  data  points  contributing  to  the 
correlation.  Under  the  assumption  that  the  two  independent 
sample  correlation  coefficients  come  from  the  same 
papulation,  their  difference  is  distributed  normally  with 
mean  0  and  variance 


2 


N  -  3 


We  can  therefore  go  to  the  standard  normal  tables  to  obtain 
the  significance  of  the  statistic 


Z1  -  Z2 


SORT  C2/(N-3)1 

which  tests  the  hypothesis  that  the  two  correlation 


coefficients  are  equal.  These  values  are  tabulated  in  Table 


r*  r* 


TABLE  4  TEST  OF  THE  HYPOTHESIS  THAT  THE 
CORRELATION  COEFFICIENTS  WITHIN  A  DATA  SET 

ARE  EQUAL 


4.  To  gat  a  significance  of  greater  than  .05  in  this  test 
requires  a  value  between  +/-  1.96,  so  only  3  of  the  18  pairs 
of  correlation  coefficients  are  statistically  equal,  giving 
very  strong  evidence  that  the  correlation  coefficients  of 
the  residuals  as  recorded  by  the  two  sensors  in  a  crossover 
pair  are,  in  general,  not  equal. 

The  practical  significance  of  the  preceding  tests  is 
that  there  appears  to  be  much  local  variation  in  the  range 
and  that  one  single  model  for  simulating  random  replications 
of  underwater  track  is  not  apparent.  Therefore,  each  case 
must  be  simulated  and  studied  separately,  and  the  study  of 
the  error  estimates  distributions  done  on  a  case  by  case 


IV.  SIMULATION  MODEL 


The  simulation  model  itself  is  a  logical  extension  of 
the  residual  generation  program.  In  -fact,  the  method  used 
to  simulate  a  data  track  was  to  simulate  residuals  and  add 
them  to  the  -fitted  straight  line  obtained  -for  that  track, 
rather  than  to  start  from  scratch  and  simulate  a  totally  new 
track  at  each  iteration.  This  method  was  very  quick  and 
yielded  very  good  simulated  track  segments.  Figure  4  on  the 
fallowing  page  is  graphical  comparison  of  a  typical  original 
track  segment  overlaid  by  a  track  segment  simulated  by  the 
method  stated  above.  The  comparison  demonstrates  the 
realistic  quality  of  the  simulated  track. 

The  regression  method  used  to  compute  residuals  of  a 
track  segment  also  produced  the  best  straight  line  in  three 
dimensional  space  to  approximate  the  target  path  through  the 
water.  Given  the  straight  line,  the  residuals  themselves, 
and  the  assumption  of  normality  of  the  residuals,  it  is 
passible  to  simulate  a  set  of  residuals  from  normal  (0,1) 
deviates  and  add  them  to  the  straight  line  to  obtain  a 
simulated  track. 

In  simulating  residuals,  the  object  is  to  compute  an 
N  x  3  matrix,  using  normal  (0,1)  deviates,  whose  vectors  of 
X,  Y,  and  Z  components  are  normally  distributed  with  mean  0 
and  whose  covariance  matrix  is  equal  to  the  covariance 


matrix  of  residuals.  That  is 

R  -  T  x  X ' 

where 

R  =  3  x  N  matrix  of  simulated  residuals 
X  =  N  x  3  matrix  of  N(0,1>  deviates 
T  =  3  x  3  transf ormation  matrix  such  that 

T  x  T'  =  covariance  matrix  of  residuals. 

Note  that  R  must  be  transposed  to  get  the  desired  N  x  3 
residual  format. 

It  is  easy  to  show  that  any  3x3  matrix  T  will  not 
alter  the  mean  of  0. 

ECR3  -  ECT  x  X'l 
=  T  x  EC X  '  3 

Since  the  expectation  of  X  —  consisting  of  normal  (0,1) 
deviates  —  is  identically  0, 

ECRU  =  T  x  (0)  =  0 
giving  the  desired  result. 

The  condition  that  T  x  T'  is  equal  to  the  covariance 
matrix  of  the  residuals  is  necessary  because 
C0VCR3  —  ECR  x  R'3/N  (the  mean  is  0) 

=  ECT  x  X'  x  X  x  T  '  3/N 
=  T  x  EC X '  x  X3/N  x  T' 

(since  T  is  a  linear  operator) 

=  T  x  I  x  T' 

(because  the  covariance  matrix  of  N(0,1)  random  variables  is 


the  identity  matrix) 


T  x  T' 


The  problem  now  becomes  one  o-f  -finding  a  3  x  3  matrix  such 
that 


COVtRl  -  T  x  T* 

We  have  the  covariance  matrix  o-f  the  residuals.  Let  T  be  an 
upper  triangular  matrix.  Then 


T11 

T12 

T13 

Til 

0 

0~ 

COVn 

COV  12 

C0V13 

0 

T22 

T23 

X 

T12 

T22 

0 

= 

cov2i 

cov22 

C0V23 

0 

0 

T33 

Il3 

T23 

Ml 

10 

1- 

f0V31 

C0V32 

C0V33^ 

Noting  that  the  covariance  matrix  is  symmetric  and 
performing  the  matrix  multiplication  yields 


T33  =  SQRT  <C0V33> 


cov 


32 


23 


33 


T22  *  SQRT  <C0V22  -  T23‘> 


COV 


13 


13 


33 


C0V21  ”  T23T13 


12 


22 


Tl,  '  SQRT(COVu  -  T132  -  Tj2Z> 


Since  the  matrix  T  is  easily  computed  and  X  can  be  generated 
from  a  random  number  generating  package,  the  residuals,  R, 
are  quickly  and  economically  computed  for  as  many  simulated 
tracks  as  desired. 


Adding  the  residuals  thus  formed  to  the  best  straight 
line  target  path  of  the  original  data  yields  a  simulated 
track. 

The  commented  FORTRAN  program  written  to  perform  the 
simulation  is  included  as  Appendix  A.  It  is  basically  a 
driver  program  that  generates  a  user  specified  number  of 
simulated  tracks  and  interfaces  with  KEYMA1N,  which 
estimates  displacement  and  angular  rotation  values.  Because 
some  of  the  user  friendly  attributes  of  KEYMAIN  interfere 
with  the  speedy  generation  of  the  simulation's  required 
output  parameters,  the  package  has  been  altered  somewhat  to 
increase  speed.  The  output  of  the  driver  program  is  two 
data  files.  One  data  file,  ROTATE.DAT,  is  an  N  x  4  matrix 
that  gives,  in  the  first  column,  the  maximum  angle  of 
rotation  of  the  sensor,  and,  in  the  last  3  columns,  provides 
the  ordered  Euler  angle  components  that  form  the  single 
maximum  rotation  angle.  The  second  file,  DISPLACE.DAT,  is 
also  an  N  x  4  matrix.  The  first  column  is  the  magnitude  of 
displacement  of  the  array,  while  the  last  3  columns  give  the 


X,  Y,  and  Z  axis  components  of  displacement. 


V.  INTERFACE  WITH  EXISTING  PROGRAMS 


The  software  developed  to  compute  residuals  and  generate 
simulated  track  segments  are  the  original  work  of  the 
author,  aided  by  such  canned  routines  as  eigensystem 
analysis,  random  number  generators  and  vector  arithmetic. 
These  canned  subroutines  came  from  the  IMSL  Library  of 
Mathematic  and  Statistical  functions  developed  for  the  IBM 
PC  computers  CRef.41,  and  are  compiled  in  machine  language 
libraries  not  reproducible  here.  The  author's  FORTRAN  code 
is  listed  as  Appendix  A. 

In  order  to  produce  the  estimates  of  sensor  displacement 
and  rotation,  however,  it  was  necessary  to  interface  with  a 
large  stand  alone  FORTRAN  package,  KEYMAIN.  This  program 
takes  as  input  crossover  region  data  sets  (real  or 
simulated)  and  produces  the  estimates  of  displacement  and 
rotation  of  the  second  sensor  of  the  crossover  pair,  based 
on  the  assumed  accurate  position  of  the  first.  The 
orientation  correction  output  by  KEYMAIN  is  actually  a  three 
valued  vector  of  ordered  angular  rotations  in  the  XY ,  XZ  and 
YZ  planes.  Similarly,  the  location  correction  is  a  three 
valued  vector  of  displacements  in  the  X,  Y  and  Z  directions. 
To  reduce  complexity,  the  ordered  Euler  angles  were  reduced 
to  a  single  maximum  angle  of  rotation  about  an  appropriately 
tilted  axis,  and  the  three  components  of  displacement  were 

40 


v.’-  v,>  >.  v  \>  .y.v 

•  ■ .v.  m  jv  wVvvi  uy  k  *  -  k  -*?■  -*%!  ^  w  - 


Ji-k.  jVi. 


reduced  to  a  single  quantity,  magnitude  of  displacement. 
Thus  the  six  dimensional  quality  o-f  position  corrections  was 
reduced  to  two. 

KEYMAIN  was  designed  as  a  stand  alone  product  and  as 
such  is  very  user  -friendly.  It  was  also  designed  to  take 
not  just  one,  but  several,  crossover  data  sets  and  produce 
displacement  and  rotation  estimates  -for  several  sensors  at 
once.  Because  its  use  in  the  simulation  was  to  process  a 
single  simulated  crossover  data  set,  and  because  it  was 
being  called  as  a  subroutine  rather  than  used  as  a  stand 
alone  package,  significant  changes  were  required  in 
the  program  to  obtain  -fast  simulated  results  uninterrupted 
by  the  now  unnecessary  user  -friendliness.  Not  only  did  this 
require  the  modi-f ication  o-f  the  executive  driver  routine, 
but  also  modification  of  several  of  the  called  subroutines. 


VI.  RESULTS 


It  is  imperative  that  the  simulated  track  segments 
"look"  and  "act"  like  the  original  track  segment  they  are 
simulating.  That  the  simulated  track  segments  act  like  the 
original  is  guaranteed  by  the  equations:  the  residuals  will 
have  mean  zero  by  definition  and  their  covariance  matrices 
were  computed  to  be  the  same  as  the  original's.  For  visual 
verification.  Appendix  D  is  included.  Appendix  D  contains 
plots  of  each  original  track  segment  overlaid  by  one  of  the 
tracks  simulated  from  it.  A  good  fit  is  apparent  in  all 
cases. 

The  simulation  model  can  produce  the  data  needed  to 
produce  a  scatterplot  of  magnitude  of  displacement  (in  feet) 
versus  maximum  angle  of  rotation  (in  radians)  like  the 
schematic  shown  in  Figure  5.  This  represents  an  imaginary 
situation  where  a  target  vehicle  was  driven  through  a 
single  crossover  region  on  a  range  700  times  over  the  exact 
same  track,  and  the  results  were  fed  into  the  program  to 
produce  the  displacement  and  rotation  values.  It  acts  as  a 
data  base,  or  sampling  distribution,  for  an  array  whose 
position  has  not  changed,  and  the  graph  depicts  the  natural 
variability  of  the  data.  The  contours  represent  some 
theoretical  confidence  levels  for  that  particular  sensor. 


For  example, 


take  the  outermost  contour,  labeled  .90 


In  a 


E  OF  DISPLACEMENT  ( 


-future  tracking  exercise,  crossover  data  can  be  -fed  into  the 
program  that  produces  displacement  and  rotation  estimates 
That  will  produce  a  point  on  the  graph.  I-f  that  point  lies 
outside  the  contour,  we  have  good  evidence  that  the  sensor 
has  moved.  Statistically  speaking,  one  would  have  a  10 
percent  chance  o-f  rejecting  a  true  null  hypothesis,  given  a 
null  hypothesis  o-f  no  sensor  movement.  Conversely,  if  the 
point  plotted  closer  to  the  middle  of  the  graph,  there  would 
be  insufficient  evidence  to  support  the  hypothesis  of  sensor 
movement,  and  the  apparent  movement  would  be  attributed  to 
the  inherent  variability  in  the  data. 

It  is  unreasonable  to  expect  that  a  target  vehicle  could 
be  driven  along  the  same  track  700  times,  nor  would  the 
range  operators  be  likely  to  attempt  it.  The  simulation 
program,  however,  needs  only  one  straight  line  segment  of 
track  through  a  crossover  region,  and  can  then  generate  as 
many  track  data  sets  as  needed  to  produce  the  graph. 

Figures  6  through  11  are  plots  produced  from  the 
simulation  model  using  the  data  sets  from  Figure  3  (p.  18) . 
Figure  9  is  a  "well-behaved”  plot  that  could  conceivably, 
with  many  more  runs,  yield  the  type  of  graph  displayed  in 
Figure  5.  In  fact,  the  points  appear  so  evenly  distributed 
that  one  could  conceivably  assume  independence  between  the 
rotation  and  displacement  values. 

Although  computationally  attractive,  independence  is 
thought  to  be  a  poor  assumption,  because  if  a  sensor  has 
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di spl aced 


a  great  distance. 


it  has  also  had  ample 


opportunity  to  rotate.  This  observation  is  speculative,  but 
seems  to  be  borne  out  in  the  remaining  figures  where  a 
strong  positive  correlation  appears  to  exist  between  the 
displacement  and  rotation  values.  The  fishhook  appearance 
on  these  graphs  is  a  result  of  making  the  correction 
estimates  positive  regardless  of  the  direction  of 
displacement  or  rotation.  The  graphs  taken  as  a  whole  serve 
to  illustrate  the  variable  nature  of  the  track  data. 

These  three  graphs,  drawn  from  the  simulation,  vividly 
illustrate  the  need  for  further  study.  Two  points 
concerning  the  data  on  which  these  graphs  were  based  are 
perhaps  of  some  importance  and  may  begin  to  explain  the 
strange  nature  of  the  graphs  produced.  First,  the  data  is 
fairly  old,  taken  in  September,  1982.  The  range  operators 
from  NUWES  state  that  their  capabilities  have  improved  in 
the  interim  3  years  and  that  newer  data  could  prove 
significantly  more  accurate.  Second,  the  data  was  taken  on 
a  single  day,  from  a  single  range,  using  only  3  of  the  24 
sensors  on  the  range. 

The  great  variety  of  graphs  produced  points  to  the 
necessity  to  do  much  more  research  and  to  the  utility  of  a 
simulation  model  to  accomplish  it.  In  addition  to  the 
aforementioned  use  of  newer  data,  the  following 
considerations  warrant  study  to  ascertain  their  possible 


effects 


Day  to  day  variations  on  tha  range  —  It  is 
not  cl aar  at  this  time  whether  the  character  of 
the  water  column  in  which  the  target  operates 
is  invariant  over  time.  Di f f erences  in  salinity 
and/or  temperature  could  conceivably  evolve  over 
time,  affecting  the  accuracy  of  the  track  data. 

Seasonal  variations  on  the  range  —  While  it 
is  not  clear  whether  changes  occur  on  the  range  on 
a  daily  basis,  it  certainly  seems  reasonable  to 
expect  a  variation  from  season  to  season.  Our 
data,  taken  from  a  single  range  on  a  single  day, 
was  insufficient  to  explore  this. 

Location  on  the  range —  Current  practice  on  a 
tracking  day  is  to  take  one  sample  of  the  water 
column  on  the  range,  checking  for  salinity, 
temperature,  and  other  factors  that  affect  the 
sound  velocity  profile,  and  assume  that  the 
results  hold  true  for  the  duration  of  the  exercise 
for  the  entirety  of  the  range.  The  possibility  of  a 
daily  change  was  discussed  earlier.  Here  we  mention 
the  possibility  of  different  water  characteri sties 
from  one  end  of  the  range  to  the  other. 

Geometry  of  tracking  runs  —  It  is  possible 
that  varying  the  geometry  of  the  tracking  run 
could  result  in  changes  to  the  quality  of  the 
data  recorded.  It  can  be  seen  from  figure  3  (on 
page  IB)  that  all  of  the  data  used  for  the 
simulations  was  from  tracks  that  run  predominantly 
downrange  with  relatively  little  crossrange  change 
and  virtually  no  depth  change.  (Although  depth 
cannot  be  seen  from  figure  3,  it  was  examined  for  all 
the  tracks. ) 

Depth  of  target  —  Depth  appears  to  be  the  least 
reliable  of  the  three  dimensions  recorded  during 
tracking  runs.  Our  data  is  from  targets  that 

operated  in  a  narrow  depth  band.  Deeper  or 

shallower  targets  could  significantly  affect  the 
data. 

Inhomogeneities  in  the  water  —  It  is  quite 

conceivable  that  undetected  inhomogeneities  in  the 
water  could  significantly  affect  the  quality  of 
the  data.  It  would  appear  that  currents  and 
turbulence  could  affect  the  passage  of  sound 

through  the  water  column,  but  quantifying  these 

disturbances  could  be  a  major  practical  problem. 


Precisely  characterizing  the  variability  of  the 
correction  estimates  is  elusive  and  the  preceding 
considerations  invite  much  future  research  before  that  goal 
is  reached.  The  existence  of  this  simulation  model 
brightens  the  outlook  for  ultimate  success. 


VII.  AREAS  FOR  FUTURE  WORK 


Problems  encountered  in  the  development  of  the 
simulation  model  and  intuition  gained  as  a  result  of  working 
on  it  gave  birth  to  several  areas  for  potential  future 
study.  Some  of  these  are  listed  below. 

A.  RESIDUALS  FOR  tURVED  CROSSOVER  DATA  TRACKS 

It  was  convenient  in  this  simulation  model  to  use 
straight  line  segments  of  track  to  get  residuals.  This  is 
because  the  method  of  principle  components  works  only  for 
straight  lines  in  N-space,  rather  than  for  curves.  If  a 
method  could  be  developed  to  regress  the  data  for  any  track 
onto  its  best  path,  regardless  of  curvature,  one  major 
obstacle  in  the  use  of  the  simulation  would  be  removed. 
Whereas  now  we  are  limited  in  the  type  of  data  we  can  use, 
such  a  method  would  enable  the  use  of  all  crossover  data. 

B.  WEIGHTED  DATA  BASED  ON  DISTANCE  TO  SENSOR  ARRAY 

It  was  assumed  in  the  model  that  recorded  track  data 
was  uniformly  accurate  and  reliable,  regardless  of  the 
distance  from  target  to  sensor.  That  is,  no  distinction  was 
made  between  the  accuracy  of  the  data  when  a  target  vehicle 
was  "close"  to  a  sensor  and  when  the  target  was  far  away. 
Some  sort  of  weighting  scheme  may  prove  beneficial  and  help 
smooth  out  the  current  data  discrepancies. 


C.  TR I SENSOR  CROSSOVER  DATA 


Although  smaller  in  size,  there  exist  areas  on  the 
range  where  a  target  may  be  tracked  simultaneously  by  three 
sensors  at  once.  These  tri sensor  crossover  regions  could 
provide  insight  into  the  accuracy  o-f  the  position  of  a 
sensor  and  may  yield  more  accurate  estimates  o-f  displacement 
and  rotation. 

D.  TRACK  DATA  SELECTION  PROCEDURES 

It  is  not  now  known  whether  a  greater  number  o-f  data 
points  in  a  track  segment  yields  better  results  in  the 
simulation  and  intuition  is  limited  on  this  point.  Research 
in  this  area  could  provide  valuable  guidelines  -for  selection 
a-f  data  to  use  in  the  simulation  in  the  -future. 


E.  DECISION  RULES  TO  DETERMINE  SENSOR  MOVEMENT 

After  accurate  character i zat i on  of  the  variability  of 
track  data  and  correction  estimates  is  made,  the  next 
logical  and  extremely  useful  step  would  be  the  construction 
of  statistically  sound  decision  rules  to  determine  the 
f ol lowing: 

(a)  The  discrepancy  noted  between  the  tracks  in  a 

crossover  data  set  can  be  explained  by  the  inherent 
variability  of  the  data. 

(b)  The  discrepancy  noted  between  the  tracks  in  a 

crossover  data  set  can  be  quantified  by  the  sensor 
slippage  model  and  can  be  computationally  corrected. 

(c)  The  discrepancy  noted  between  the  tracks  in  a 

crossover  data  set  cannot  be  explained  by  the  model 
of  sensor  movement  and  some  other  explanation  must  be 
sought. 


Clearly,  the  last  option  is  the  least  desirable,  but  if 
that  situation  is  present,  it  needs  to  be  noted. 

F.  COMPUTATIONAL  RANGE  RESURVEY 

If  the  method  of  providing  correction  estimates  can  be 
proven  to  be  reliable  and  accurate,  it  provides  a  potential 
method  of  range  resurvey  that  has  several  advantages  over 
the  current  method.  First,  the  range  would  not  need  to  be 
shut  down  to  resurvey.  In  fact,  range  use  would  be 
mandatory  to  keep  up  to  date  sensor  array  positions. 
Second,  it  would  be  less  expensive  than  the  current  method, 
replacing  the  equipment  and  manpower  intensive  current 
process  with  relatively  inexpensive  computer  assets.  Third, 
it  would  take  less  time.  Resurvey  of  a  single  sensor  array 
on  the  range  can  take  up  to  a  day;  generating  corrections 
from  track  data  takes  seconds.  Fourth,  since  the  current 
method  uses  a  craft  on  the  surface  equipped  with  a  pinger, 
all  the  pings  must  travel  through  the  first  150-200  feet  of 
the  water  column,  where  the  sound  velocity  profile  is  quite 
variable  and  most  difficult  to  determine,  to  get  to  the 
sensor  array.  In  contrast,  the  underwater  target  vehicles 
tracked  by  the  arrays  are  typically  in  400-600  feet  of  water 
where  the  sound  velocity  profile  is  much  smoother  and  easier 
to  predict.  Thus,  one  source  of  variability  in  determining 
sensor  position  is  reduced. 

It  is  not  envisioned  that  this  computer  method  could 
ever  replace  the  current  survey  process,  but  rather  augment 


it.  Soma  way  to  determine  tha  position  o-f  ona  sansor  is 
necessary  before  KEYMAIN  can  even  begin  to  function. 
However,  if  this  computer  process  could  augment  current 
resurvey  efforts,  or  reduce  the  frequency  with  which 
resurveys  must  be  made,  significant  savings  in  time  and 
money  could  result. 


APPENDIX  A 


tiiV : 


FORTRAN  LISTING  FOR  PROGRAM  SIMDAT 


PROGRAM  SIMDAT2 

This  program  simulates  3-D  track  data  based  on  a 
specific  real  track  specified  by  the  user. 

User  inputs  are: 

1.  track  segment  data  file  to  be  simulated 

2.  number  of  the  left  and  right  sensor  in  the 
crossover  pair 

3.  number  of  simulated  tracks  desired 

4.  request  for  sample  simulated  track  (YES  or 

5.  random  number  generator  seed 


The  program  output 


file  of  residuals  from  the  original  track 
( RES I DUAL . DAT ) 

file  of  a  simulated  crossover  track,  if 
requested  (SIMTRACK.DAT) 

file  of  displacement  values  (in  feet)  for  the 
right  sensor  of  each  simulated  track  in  4 
columns 

Col  1  magnitude  of  displacement 
Col  2-4  X,Y,Z  components  of  displacement 
file  of  rotation  values  (in  radians)  for  the 
right  sensor  of  each  simulated  track  in  4 
col umns 

Col  1  maximum  angle  of  rotation 

Col  2-4  ordered  Euler  angles  of  rotation 


SIMS, 


VARIABLE  DECLARATION 

INTEGER*4  N,  I,  J,  K,  IER,  BIG1 (3) ,  BIG2(3) 
+  POINT ( 130) ,  TRACKS,  IDL,  IDR,  TRKOUT 

CHARACTER  DSNAME*13 

REAL*4  NORM (260) 


REAL*8  TRACK ( 130,6) ,  TR1 (130,3),  TR2<130,3),  MBAR 1(3), 
i-  MBAR2  (3)  ,  MBARM1  (130,3)  ,  MBARM2  ( 130, 3)  ,  C0V1<3,3), 

-  C0V2 (3,3) ,  DIF1 ,  PI (3,3),  MUT1 (130,3),  S1M1  CS, 1 30)  , 

h  DD1  (3)  ,  DD2  (3)  ,  PA(3,3),  PB<3,3),  W0RM130),  L>1  ,  D2, 
■  DIF2,  SUM1 ,  SUM2 ,  A1 ,  A2,  Zl(3,130>,  Z2(3,130>, 
h  P2 (3,3) ,  ZT1 (130,3),  ZT2<130,3),  TBAR ,  Z1BAR,  Z28AR, 

-  f RKSUM ( 130) ,  TKSUM2,  CT1(3,130),  CT2(3,13U),  SEED, 

-  MU  T2  < 1 30 , 3)  ,  RESID1 ( 130,3)  ,  RESID2 ( 130,3)  ,  SU2<3,3), 
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+  CO V 1 R  <  3 , 3 ) ,  C0V2R<3,3),  SQ1<3,3>,  SIMTRK <200,6) , 

+  SI M2 (3, 130) ,  ROTATE < 1000,4) ,  DISP < 1000,4) ,  DAIA<2,4> 

...Begin  the  user  input  section 

WR I TE  ( * , * )  'Enter  FNAME.FT  of  crossover  data  set 
WRITE <*. #)  '  on  disk  s  ' 

READ  <*, ' <A)  *  >  DSNAME 
WRITE <*,*)  '  ' 

WRITE  <*,*)  'What  is  the  NUMBER  o-f  the  left  sensor  in 
WRITE <*.*)  '  the  crossover  pair  ?  ' 

READ  <*,*)  IDL 
WRITE  <*,*)  '  ' 

WRITE<*,*)  'What  is  the  NUMBER  of  the  right  sensor  ? 
READ  <*, *)  IDR 
WR I TE  <  * , * )  '  ' 

WRITE (*,*>  ‘How  many  simulated  tracks  do  you  desire  ? ' 
WRITE <*.*>  '(NOTE  s  max  1000)  ' 

READ(*,*>  TRACKS 
WRITE  <*,*)  '  ' 

WRITE  <*,*)  'Do  you  want  a  sample  simulated  track  ?' 
WRITE  <  * , *)  'Enter  1  for  YES,  0  <zero>  for  NO 
READ  <  * , * )  TRKOUT 
WRITE  <*,  *)  '  ' 

WRITE<*,*)  'Seed  for  the  random  number  generator’ 

WRITE <*,*)  'NOTE  :  Seed  must  include  a  decimal 
READ (*, *)  SEED 
WRITE  < * , #)  '  ' 

...Read  data  from  file 


OPEN  <  1  ,  F I  LE=*DSNAME ,  STATUS*  '  OLD  '  ) 

N  =  O 

10  N  =  N  +  1 

READ ( 1 , * , END=30 ,  ERR=30>  F OINT <N) , (TRACK <N , I ) , I = 1 , 6) 

...Separate  into  "left"  and  "right"  sensor  tracks 

DO  20  I  *  1,3 

TR1 <N, I)  =  TRACK  <N , I ) 

TR2  <N , I )  =  TRACK  <N, 1+3) 

20  CONTINUE 
GOTO  10 

30  CLOSE  (UNIT  =  1) 

N  =  N  -  1 
C 

C. . . Compute  the  covariance  matrix  for  each  track. 

C... -First  step,  get  column  averages  (with  first  column 
C. . .  average,  TBAR,  computed  for  later  use) 

C 

TBAR  =  0. 

DO  50  I  =  1,3 


-  ^  ■ 
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MBAR1 (I)  >0. 

MBAR2 < I )  *  0. 

DO  40  J  =■  1,N 

MBARKI)  ■  MBARl(I)  +  TR1(J,I) 

MBAR2  < I )  »  MBAR2  < I )  +  TR2(J,I) 

IF  (I  .EQ.  1)  TBAR  =  TBAR  +  DBLE (P0IN1 < J ) ) 

40  CONTINUE 

MBARKI)  =  MBARKI)  /  DBLE  (N) 

MBAR2  < I )  =  MBAR2  < I )  /  DBLE (N) 

50  CONTINUE 

TBAR  ■  TBAR  /  DBLE (N) 

C 

C...Do  the  matrix  multiplication  t  A(t)xA  ,  subtracting  off 
C...the  mean  from  each  column  entry  to  form  the  covariance 
C. . . matrix.  Also  make  the  matrix  of  (points  -  means)  for 
C.  . . 1 ater  use. 

C 

DO  80  I  =  1,3 
DO  70  J  ■  1,3 
C0V1 <I,J>  =  0. 

CQV2 < I , J )  «  0. 

DO  60  K  =  1  ,  N 

C0V1 < I , J )  =  COV 1(  I , J )  +  ( TR 1( K , I ) -MBAR 1(1)) 

*  * ( TR 1 ( K , J ) -MBAR 1 ( J ) ) 

C0V2 ( I , J )  =  CO  V2 ( I , J ) + ( TR2 ( K ,  I ) — MBAR2 ( I ) ) 

*  * ( TR2  <  K , J ) -MBAR2 ( J ) ) 

60  CONTINUE 

C0V1 ( I , J )  =  COV 1 ( I , J )  /  DBLE(N-l) 

C0V2 ( I ,  J )  =  C0V2 ( I ,  J )  /  DBLE(N-l) 

70  CONTINUE 
80  CONTINUE 

DO  lOO  I  »  1,N 
DO  90  J  =  1,3 

MB ARM 1 ( I , J >  *  TR1(I,J)  -  MBARl(J) 

MBARM2 ( I , J )  =  TR2 ( I ,  J  )  -  MBAR2 ( J ) 

90  CONTINUE 

100  CONTINUE 
C 

C. . .Form  the  matrix  P  of  ordered  principle  components  for 
C...each  track.  The  columns  of  P  are  the  eigenvectors 
C. .. assoc i at ed  with  the  eigen-values  of  the  covariance 
C... matrix  for  each  track  arranged  in  order  of  descending 
C. . .eigenvalues.  (ie.  the  eigenvector  associated  with 
C...the  largest  eigenvalue  is  the  first  column) 

C 

C...Call  routine  to  compute  eigenvalues/vectors  for  each 
C. . . track 
C 

CALL  E I GRS ( COV 1,3,1 1 ,DD1 , PA ,3, WORK , I ER ) 

CALL  E I GRS  <  C0V2 ,3,1 1 , DD2 , PB , 3 , WORK , IER) 

...Get  eigenvectors  in  eigenvalue  order,  largest  to 
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smallest,  by  column 


DD1/2  =  eigenvalue  vector 
PA/PB  »  eigenvector  matrices 

BlGl(l)  -  1 
BIG2 ( 1 )  -  1 
BIG1 (3)  -  3 
BIG2 <3)  -  3 

IF  (DD1 (BIG1 ( 1 ) )  .LT.  DD1 (2) )  BIGl(l)  =  2 

IF  <DD2<BIG2<D )  .LT.  DD2<2>>  BIG2<1)  =  2 

IF  <DD1 (BIGl <1) )  .LT.  DD1 <3) )  BIG1<1>  =  3 

IF  <DD2<BIG2<1D  .LT.  DD2<3>)  BIG2<1)  -  3 

IF  <DD1<BIG1<3D  ,  GT .  DDK1D  BIGl  <3)  -  1 

IF  <DD2<BIG2<31 )  .GT.  DD2 < 1 D  BIG2<3>  =  1 

IF  <DD1<BIG1<3D  .GT.  DD1  (2D  BIGl  (3)  -  2 

IF  <DD2<BIG2<3D  .GT.  DD2<2D  BIG2<3>  *  2 

IF  <  <BIG1 < 1 )  +  BIGl <3) )  .EQ.  3)  THEN 
BIGl <2)  «  3 
ELSE 

IF  <  (BIGl  <  1 )  +  BIG1<3D  .EQ.  4)  THEN 
BIGl <2)  *  2 
ELSE 

BIGl <2)  =  1 
END  IF 
END  IF 

IF  < (BIG2 < 1 )  +  BIG2  <3) >  .EQ.  3)  THEN 
BIG2  <2)  -  3 
ELSE 

IF  <  < B I G2  < 1 )  +  B I G2  < 3 ) )  .EQ.  4)  THEN 
BIG2  <2)  =  2 
ELSE 

BIG2  <2)  -  1 
END  IF 
END  IF 

DO  130  1-1,3 
DO  120  J  -  1,3 

PI <1  ,J)  *  PA  < I , BIGl <  J  D 
P2  < I , J )  *  PB  < I , BIG2 ( J ) ) 
i  CONTINUE 

i  CONTINUE 

Compute  the  matrix  ZT  for  each  track 

ZT  represents  the  projection  of  the  track  data  onto 

the  principle  components 

ZT  =  <TR  -  MBAR)  x  P  -  MBARM  x  P 

where  TR  -  MBAR  =  track  data  minus  the  column  average 
for  each  row 

Call  routine  to  multiply  matrices  AxB 


CALL  VMULFF  < MBARM 1 ,P1 , N , 3 , 3 , 130 , 3 , ZT1 , 130, IER) 


CALL  VMULFF ( MBARM2 , P2 , N , 3 , 3  f 1 30 , 3 , ZT2 , 1 30  , I ER ) 


...Since  there  are  some  points  missing  -from  the  data  set, 
.  ..perform  a  simple  least  squares  linear  regression  onto 
...the  first  principle  component 
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160 


TKSUM2  =0.0 
Z1BAR  =  0.0 
Z2BAR  =0.0 
DO  140  I  »  1,N 

Z1BAR  =  Z1BAR  +  ZT1(I,1) 

Z2BAR  =  Z2BAR  +  ZT2(II1) 

TRKSUM  < I )  =  DBLE (POINT ( I >  >  -  TBAR 
TKSUM2  =  TKSUM2  +  TRKSUM(I)#*2 
CONTINUE 

Z1BAR  =  Z1BAR  /  DBLE (N) 

Z2BAR  =  Z2BAR  /  DBLE (N) 

SUlil  =0.0 
SUM2  =0.0 
DO  150  I  =  1,N 

DIF1  =  (ZT1 (1,1)  -  Z1BAR)  *  TRKSUM(I) 

DIF2  =  ( ZT2 (1,1)  -  Z2BAR)  *  TRKSUM < I ) 

SUM1  =  SUM1  +  DIF1 

SUM2  =  SUM2  +  DIF2 

CONTINUE 

D1  =  SUM1  /  TKSUM2 

D2  =  SUM2  /  TKSUM2 

A1  =  Z1BAR  -  D1  *  TBAR 
A2  =  Z2BAR  -  D2  *  TBAR 
DO  160  I  =  1,N 

ZT1 (1,1)  =  A1  +  D1  *  DBLE (POINT ( I ) ) 
ZT2  (1,1)  =  A2  +  D2  *  DBLE (POINT  < I ) ) 
CONTINUE 


...Get  CT  matrix  which  represents  the  orthogonal  projection 
...of  the  data  onto  the  straight  line  of  the  first 
.  . .principle  component 
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ISO 


DO  170  I  = 
Z1 (1 , I) 
Z2 (1,1) 
CONTINUE 
DO  180  I  = 
Z1 (2,1) 
Z1 (3,1) 
Z2(2, I) 
Z2 (3 , I ) 
CONTINUE 


1  ,  N 

=  ZT1  (1,1) 
=  ZT2 (1,1) 


1,N 
=  0.0 
=  0.0 
=  0.0 
=  0.0 


...Call  routine  to  multiply  matrices  AxB 


CALL  VMULFF (P1,Z1,3,3,N,3,3, CT 1 ,3, IER) 


CALL  VMULFF  (P2 , Z2 , 3 , 3 , N , 3 , 3 , CT2 , 3 , IER) 


...Move  "line"  of  data  back  into  original  coordinate  system 

DO  200  I  -  1,N 
DO  190  J  -  1,3 

MUT 1 ( I , J )  -  CT 1  ( J  ,  I )  +  MBAR1 < J ) 

MUT2  < I , J )  -  CT2  < J , I )  +  MBAR2 ( J ) 

...Compute  residuals  for  each  track 

RES1D1 ( I , J)  -  TR1 ( I ,  J )  -  MUT 1  ( I ,  J ) 

RESID2 ( I , J )  -  TR2  <  I ,  J)  -  MUT2(I,J) 

190  CONTINUE 

200  CONTINUE 

...Write  the  set  of  residuals  out  to  the  file  RESIDUAL.DAT 


OPEN  <2,  FILE  -  'RESIDUAL.DAT',  STATUS  =  'NEW') 

DO  210  I  ■  1,N 

WRITE (2,330)  (RESIDl (I ,J) , J=1 ,3) , (RESID2 ( I , J) ,J=1 ,3) 
210  CONTINUE 

...Compute  the  covariance  matrix  of  the  residuals 
...(Note  :  column  averages  are  identically  zero) 

...Call  routine  to  multiply  matrices  A(t)xB,  then  divide 
.  . . by  N-l 

CALL  VMULFM ( RES I D 1 , RESIDl ,N,3,3, 130, 130,C0V1R,3, IER) 
CALL  VMULFM ( RES I D2 , RES I D2 , N , 3 , 3 , 1 30 , 130,C0V2R,3, IER) 

DO  230  1=1,3 
DO  220  J  =  1,3 

C0V1R ( I , J )  =  COV 1 R ( I , J )  /  DBLE(N-l) 

C0V2R ( I , J )  =  C0V2R ( I , J >  /  DBLE(N-l) 

220  CONTINUE 

230  CONTINUE 


.Get  "square  root"  of  covariance  matrix  of  residuals 


SQ1 

(3 

,3) 

SQ2 

(3 

,3) 

SQ1 

(2 

,3) 

SQ2 

<2 

,3) 

SQ1 

(2 

,2) 

SQ2 

(2 

,2) 

SQ1 

(1 

,3) 

SQ2 

( 1 

,3) 

SOI 

( 1 

,2) 

SQ2 

( 1 

,2) 

SQ1 

( 1 

,1) 

SQ2 

( 1 

,1) 

DSQRT ( COV 1 R  <  3 , 3 ) ) 

DSQRT  <  C0V2R (3,3) ) 

C0V1R (2,3)  /  SQ1 (3,3) 

C0V2R (2,3)  /  SQ2 (3,3) 

DSQRT ( COV 1 R ( 2 , 2 )  -  SQ1 (2,3) **2> 

DSQRT ( C0V2R (2,2)  -  SQ2(2,3>**2' 

COV 1 R ( 1 , 3 )  /  SQ1 (3,3) 

C0V2R (1,3)  /  SQ2 (3,3) 

( COV 1 R ( 1 , 2 )  -  SQ1 (1,3>*SQ1 (2,3) >  /  SQ1(2,2> 
(C0V2R (1,2)  -  SQ2 (1,3) *SQ2 (2,3) )  /  SQ2(2,2) 
DSQRT (C0V1R ( 1 , 1)-(SQ1 (1 ,2>**2+SQl (1,3) **2> ) 
DSQRT (C0V2R (1,1) - (SQ2 (1,2) **2+SQ2 (1,3) **2> ) 


.  V 'V V.V 
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SQl <2, 1)  -  0. 

SQ2 (2,1)  -  0. 

SQl <3, 2)  =  0. 

SQ2 (3,2)  -  0. 

SQl (3,1)  *  0. 

SQ2 (3,1)  =  0. 

...Compute  sets  of  residuals  and  get  rotation/displacement 
. . . values 

DO  410  SIMS  -  1, TRACKS 

...Compute  set  of  simulated  residuals  from  normal (0,1) 

.  . . devi ates 

DO  260  1-1,3 

...Call  routine  to  generate  Normal (0,1)  deviates 

CALL  GGNPM( SEED, 2*N, NORM) 

DO  250  J  =  1,N 

RESID1 (J, I)  -  NORM ( J ) 

RESID2 ( J , I )  *  NORM (N+J  > 

250  CONTINUE 

260  CONTINUE 

...Call  routine  to  multiply  matrices  AxB(t) 

CALL  VMULFP (SQl , RES ID1,3,3,N,3,130,SIM1,3,I ER ) 

CALL  VMULFP ( SQ2 , RES I D2 ,3,3,N,3,1 30 , S I M2 ,3,1 ER ) 

...Put  together  Nx6  matrix  of  simulated  tracks  for  both 
.  ..arrays  by  adding  straight  line  in  original  coordinate 
...system  to  residuals 

DO  280  I  =  1 , N 

DO  270  J  -  1,3 

SIMTRK ( I , J )  =  SIM1  ( J  ,  I )  +  MUT 1 ( I , J ) 

SIMTRK ( I , J+3)  =  SI M2  ( J  ,  I )  +  MUT2(I,J) 

270  CONTINUE 

280  CONTINUE 


...Write  the  first  simulated  track  out  to  the  file 
...SIMTRACK.DAT  if  a  sample  simulated  track  was  requested. 

IF ((SIMS  .EQ.  1)  .AND.  (TRKOUT  .GT.  0) )  THEN 

OPEN  (3,  FILE  =  'SIMTRACK.DAT',  STATUS  =  'NEW') 

DO  285  I  =  1  ,  N 

WRITE  (3,340)  (SIMTRK  ( I  ,<! )  ,  J  =  1,6) 

285  CONTINUE 

END  IF 
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C...Feed  the  simulated  track  into  KEYMAIN  to  get  rotation 
C...and  displacement  numbers 
C 

CALL  KEYSUB ( S I MTRK , N , D AT A , I DL , I DR ) 

...Make  2  matrices  -  one  -for  displacement  data  and  one  -for 
...the  rotation  data 

DO  290  1-1,4 

DISP (SIMS, I )  -  DATA  (1,1) 

ROTATE (SIMS, I)  »  DATA (2, I ) 

290  CONTINUE 

WRITE  (* , 320)  SIMS 

...Go  back  and  do  it  again 

410  CONTINUE 

...A-fter  TRACKS  simulated  tracks,  write  the  displacement 
...sets  and  the  rotation  sets  out  to  a  file 

OPEN  (4 , FILE  -  'DISPLACE.DAT' , STATUS  *  'NEW') 

OPEN <5, FILE  *  'ROTATE.DAT' , STATUS  =  'NEW') 

DO  300  I  *  1, TRACKS 

WRITE (4,310)  (DISP ( I , J) , J  =  1,4) 

WRITE  (5,310)  (ROTATEd  ,  J)  ,J  *  1,4) 

300  CONTINUE 


...Close  out  the  files 

CLOSE (UNIT  -  2) 

CLOSE (UNIT  -  3) 

CLOSE (UNIT  -  4) 

CLOSE (UNIT  -  5) 

STOP 

310  FORMAT (2X,4F17. 8) 

320  FORMAT (2X, 'Through  KEYMAIN  ',16,'  time(s)  so  far’) 
330  FORMAT (IX, 6F 12.7) 

340  FORMAT (IX, 6F 12.5) 

END 
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APPENDIX  E 


FORTRAN  SUBROUTINES  CONECT  AND  REDUCE 

The  program  KEYMAIN  requires  that  any  array  -for  which 
correction  estimates  are  desired  must  be  "connected"  to  the 
first  input  sensor  array  by  one,  or  a  series  of,  crossover 
data  sets.  For  example,  if  an  input  crossover  data  set  uses 
sensor  arrays  5  and  6,  while  another  uses  arrays  6  and  10, 
all  three  of  the  arrays  (5,  6  and  10)  are  connected.  If, 
however,  a  first  input  data  set  uses  arrays  7  and  9,  a 
second  input  data  set  uses  12  and  10,  while  a  third  input 
data  set  uses  9  and  13,  the  arrays  7,  9  and  13  are 
connected,  12  and  10  are  connected,  but  all  five  arrays  do 
not  form  a  single  connected  set.  In  this  case,  if  7  was  the 
first  array  input  to  the  program,  correction  estimates  could 
not  be  made  for  arrays  10  and  12.  The  subroutine  C0NECI 
checks  to  see  that  connectedness  exists  in  the  input  data 
before  KEYMAIN  is  allowed  to  continue. 

KEYMAIN  allows  3  options  if  CONECT  discovers  that  the 
arrays  of  the  input  data  sets  are  not  connected.  One  option 
is  to  quit,  in  which  case  the  program  terminates.  Another 
option  is  to  add  more  data  so  that  all  the  arrays  are 
connected.  In  the  second  example  above,  for  instance,  if  a 
crossover  data  set  using  arrays  9  and  12  was  input  into  the 
program,  all  5  arrays  would  then  be  connected  and  KEYMAIN 


96 


Mould  continue.  A  third  option  is  to  continue  with  KEYMAIN, 
but  use  only  the  first  connected  set,  that  is,  all  the 
arrays  that  are  connected  to  the  left  array  of  the  first 
input  crossover  data  set.  This  option  presents  special 
problems  because  it  requires  a  reduction  of  the  data 
structures  that  have  been  built  as  each  data  set  was  input. 
The  subroutine  REDUCE  does  this  data  structure  reduction  and 
is  called  from  KEYMAIN  only  when  this  third  option  is 
selected. 

KEYMAIN  passes  to  CONECT  two  pieces  of  information. 
The  first  is  the  variable  Ri  that  represents  the  number  of 
data  sets  that  were  input  into  KEYMAIN.  The  second  piece  of 
information  is  a  2  x  Ri  matrix,  IND2,  that  contains  the 
number  of  the  left  and  right  sensors  of  the  RI  crossover 
data  sets.  Row  1  contains  the  left  sensor  numbers,  row  2 
the  right.  CONECT  performs  its  connectedness  check  by 
starting  a  variable  length  list  that  contains  the  array 
numbers  of  those  arrays  that  are  connected.  The  list  starts 
with  only  two  entries,  the  left  and  right  sensor  arrays  of 
the  first  crossover  data  set.  These  are  connected,  and  the 
left  sensor  is  the  "root"  to  which  all  arrays  should 
connect.  Elements  are  added  to  the  list  by  sequencing 
through  the  list  from  the  beginning,  and  adding  to  the  list 
any  array  numbers  for  IND2  that  <1)  are  not  yet  on  the  list 
and  (2)  are  connected  by  a  crossover  data  set  to  an  array 
that  is  already  on  the  list.  If,  after  sequencing  through 
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the  variable  length  list,  there  are  arrays  remaining  in  IND2 
that  never  were  put  on  the  list,  those  arrays  were  not  a 
part  of  the  first  connected  set.  If  all  arrays  in  IND2  were 
included  on  the  list,  all  the  arrays  were  connected.  After 
the  first  list  is  exhausted,  CONECT  repeats  this  procedure 
as  many  times  as  there  are  disjoint  sets  of  arrays,  starting 
each  new  list  with  the  arrays  of  a  data  set  not  in  any 
previous  set.  CONECT  informs  the  user  of  the  individual 
sets  of  connected  sets  of  arrays  and  raises  a  flag  to  alert 
the  user  if  all  sets  were  not  connected. 

If  all  the  arrays  in  the  input  crossover  data  sets  were 
not  connected,  the  user  has  three  options  to  proceed, 
described  earlier.  If  he  chooses  to  continue  using  the 
first  set  of  connected  arrays,  REDUCE  is  called  to  pare  the 
data  strauctures  built  up  during  the  data  input  process.  In 
particular,  the  R1  x  3  x  3  array  CROSSA  and  the  R1  x  6 
matrix  mean  need  to  be  reduced  to  contain  only  those 
elements  that  correspond  to  the  data  sets  in  the  first 
connected  set.  CROSSA,  which  has  R1  "pages"  of  3  x  3 
matrices  of  crossproduct  deviations  from  the  mean,  needs  to 
have  those  pages  removed  that  correspond  to  every  crossover 
data  set  in  the  original  input  not  connected  to  the  first 
array.  MEAN,  which  has  R1  rows  of  the  colunm  averages  for 
each  data  set,  needs  to  have  those  rows  removed  that 
correspond  to  data  sets  not  connected  to  the  first  array. 


The  variable  R1  itself  will  ba  raducad  to  reflect  this 


smaller  subset  of  connected  array  pairs. 

CCNECT  stores  the  information  that  indicates  which  data 
sets  in  IND2  are  connected  to  the  first  array.  This 
information  is  passed  to  REDUCE  through  KEYMAIN.  REDUCE 
reforms  the  data  structures  to  reflect  the  smaller  number  of 
data  sets  now  being  considered  by  KEYMAIN. 

CONECT  also  provides  the  matrix  IND1,  a  K  x  R1  matrix 
that  is  used  elsewhere  in  KEYMAIN.  The  R1  columns  represent 
the  crossover  data  sets  input  into  KEYMAIN.  The  variable  K 
represents  the  number  of  individual  arrays  in  the  R1  data 
sets,  each  row  representing  a  separate  array.  For  each 
column  in  IND1,  the  entries  are  all  0  except  in  the  rows 
that  represent  the  left  and  right  sensor  for  that  column's 
data  set.  If  the  row  corresponds  to  the  left  array,  the 
value  is  1.  If  it  corresponds  to  the  right  array,  the  value 
is  2.  If  REDUCE  is  called,  some  columns  (representing 
crossover  data  sets)  and  rows  (representi ng  individual 
sensors)  of  IND1  need  to  be  omitted.  REDUCE  reforms  IND1  to 
its  smaller  size. 


APPENDIX  F 

FORTRAN  LISTING  FOR  SUBROUTINE  CONECT 


SUBROUTINE  CONECT  (OUT ,R1  , IND2,K, IND1 , IA,TESTC, IND2R, 
*  DATSET) 

C 

C***************************»******************************* 
C  This  subroutine  checks  -for  the  connectedness  of  the 
C  input  data  sets.  If  the  problem  is  connected  then  the 
C  user  is  informed  and  the  array  pairs  are  printed  on  the 
C  screen;  if  not  connected,  then  the  user  is  prompted  to 

C  select  one  of  three  options  -  quit,  add  conecting  data 

C  sets,  or  run  the  program  using  the  first  connected  set 
C  that  was  input.  Gygax  -  July  1985 

C*********************************************************** 

c 

C  ...Variable  declarations. 

C 

I NTEGER*4  R1 ,K, IND2 (2,30) ,IND1 (30,30) ,I,J,IA(30> ,FIRS1 
INTEGER*4  LIST(30> , BEGIN, HALT , DISCON, L,M,0, TESTC,0U I , 

I NTEGER*4  DATSET ( 30  > , COUNT , SAVE (2,30),  I ND2R (2,30) 

C 

C  ...Initialize  the  values  of  FIRST  and  COUNT: 

C 

FIRST  =  0 
COUNT  =  0 
C 

C  ...Make  vector  IA  »  list  of  all  arrays  (w/o  repeats) 

C  in  IND2  and  get  the  value  for  K  »  #  of  individual 

C  arrays. 

C 

I A ( 1 )  -  IND2 (1,1) 

IA (2)  =  I ND2 (2,1) 

K  -  3 

IF  (R1  .EQ.  1)  GOTO  60 
DO  50  I  -  1 ,R1 
DO  40  J  ■  1,2 
M  =  K  -  1 

DO  30  L  ■  1 ,1*1 

IF  ( IND2 ( J , I )  .EQ.  I A (L) )  GOTO  40 
30  CONTINUE 

IA (K)  =  IND2 ( J , I ) 

K  =  K  +  1 

40  CONTINUE 

50  CONTINUE 
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WRITE < OUT,*) 'Rl  ' ,R1 
WRITE (OUT,*) 'K  ' ,K 

...For  each  column  of  IND1  (columns  correspond  to 
data  sets)  the  entries  are  all  zero  except  tor 
the  row  that  corresponds  to  the  left  array  1) 
and  the  right  array  <=  2) . 

DO  80  I  =  1 , Rl 
DO  70  J  =  1  ,K 
IND1  (J, I)  =0 

IF  ( IND2 (1,1)  .EQ.  I A ( J ) )  IND1(J,I>  =  1 

IF  ( I ND2 (2,1)  .EQ.  IA(J>)  IND1(J,I)  =  2 

70  CONTINUE 

30  CONTINUE 

...Check  to  see  if  all  the  arrays  are  connected. 

TESTC  =  1 

LIST ( 1 )  =  — I A ( 1 ) 

DO  131  I  =  1,R1 

IF  ( I ND2 (1,1)  .EQ.  —LIST ( 1 ) )  IND2(1,I>  =  -IND2 ( 1 , 

IF  ( IND2  (2,1)  .EQ.  -LIST(l))  IND2(2,I)  -=  -IMD2(2, 

131  CONTINUE 
BEGIN  =  1 
HALT  =  1 

140  IF  (.NOT.  (BEGIN  .LE.  HALT))  GOTO  170 
NODE  =  LIST (BEGIN) 

BEGIN  =  BEGIN  +  1 
DO  150  I  *  1 ,R1 

IF  ( .  NOT.  (  (NODE.  EQ.  IND2(1  ,  I  )  >  .  AND.  ( IND2  (2,  I  >  .  (?)  .  >  •  ' 
*  GOTO  150 

HALT  =  HALT  +  l 
LIST  (HALT)  —  IND2  (2,1) 

DO  141  J  =  1 , R 1 

IF  (I  ND2  (  1  ,  J  )  .  EQ .  —  L  1ST  ( HALT )  )  I ND2  (  1  ,  J  >  ■=-  l  MD2  (  t  .  ’ 

IF  (  IND2  ( 2 ,  J  )  .EQ. -LIST  (HALT)  )  I ND2  ( 2  ,  .7  ?  -  I  NO 

1 4 1  CUNT  I  HUE 
550  CONTINUE 

DO  1 60  I  =  1 , R 1 

IF  (  .  NOT.  (  (NODE.  EQ.  IMD2  (2,  (  >  >  .  AND.  v  IND2  \  1  ,  I  *  .  i  .  ' 

-*  GOTO  160 

HALT  HALT  +  1 
LIST '  HALT)  =  - IND2  < 1,1) 

00  131  J  =  1 , R 1 

IF  (I  ND2  (  1  ,  J )  .  EQ .  -L I  ST  ( HAL  T  >  )  I ND2  (  1  ,  J  •  ---  l  NO'..’ '  ! 

IF  ( IND2  (2,  J  )  .  EQ.  -LIST  (HALT)  )  IND2  (2  ,  J  )  -- I ND2  '  .1  . .. 

151  CONTINUE 
5  CONTINUE 
GOTO  140 
t’«!  CONTINUE 


ft  N 


DISCON  =  0 
WRITE (OUT ,230) 

DO  200  I  =  1 , R1 

IF  ( IND2 (1,1)  .LT.  0)  1T0  190 

IF  < IND2 (1,1)  .EQ.  0)  Gl  TO  200 

IF  ( < IND2 (1,1)  .GT.  0)  .AND.  (DISCON  .EQ.  1 ) ) 

FIRST  =  FIRST  +  1 

DISCON  =  1 

TESTC  =  0 

BEGIN  =  1 

HALT  =  1 

LIST ( 1 )  =  IND2 (1,1) 

GOTO  200 

190  WR I TE ( OUT , 240 )  - I ND2 (1,1) , - I ND2 (2,1) 

IF  (  (FIRST. EQ.O) .OR.  ( (FIRST. EQ. 1) .AND.  (DISCON. 
*  THEN 

COUNT  =  COUNT  +  1 
I ND2R ( 1 , COUNT )  =  -IND2(1,I> 

I ND2R ( 2 , COUNT )  =  -IND2(2,I> 

DATSET (COUNT)  =  I 
END  IF 

SAVE (1,1)  =  — IND2 (1,1) 

SAVE (2,1)  =  -IND2 (2,1) 

I ND2 (1,1)  =  O 
IND2 (2, I )  =  O 
200  CONTINUE 

IF  (DISCON  .EQ.  1)  GOTO  140 
DO  220  I  =  1 , R 1 

IND2 (1,1)  =  SAVE (1,1) 

IMD2 (2, I )  =  SAVE (2, I ) 

220  CONT I NUE 
RETURN 

30  FORMAT ( IX , 'THE  FOLLOWING  PAIRS  ARE  CONNECTED  : 
40  FORMAT ( 1 X , 1415) 

END 


GOTO  200 


EQ.  1  )  ) 
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APPENDIX  G 


FORTRAN  LISTING  FOR  SUBROUTINE  REDUCE 


SUBROUTINE  REDUCE  ( CROSS A, MEAN, Rl ,K,IND1,IA, 

*  I ND2R  ,  DATSET ) 

C 

C*********************************************************** 
C  This  is  a  specialized  subroutine  that  is  used  when 
C  option  three  is  invoked  as  a  result  of  a  failed  con- 
C  nectedness  test.  The  disconnected  data  sets  must  be 

C  removed  from  the  variables  CROSSA  and  MEAN,  and  other 

C  program  supporting  variables  must  be  adjusted. 

C  Gygax  -  July  19B5 

C*********************************************************** 
...  Variable  declarations. 

INTEGER*4  R1,K,IND1 (30,30) ,IA(30> ,IND2R(2,30> ,I,J,L, 

2  M, DATSET (30) 

REAL*8  CROSSA (30,3,3), MEAN (30,6) 

...  Compute  the  new,  reduced  Rl: 

DO  10  I  *  1,30 

IF  ( IND2R (1,1)  .EQ.  0)  GOTO  20 
10  CONTINUE 
20  Rl  =  I  -  1 

...  Make  new,  reduced  vector  I A  ■  list  of  all  arrays 
in  IND2R  w/o  repeats.  Also,  compute  a  new  K. 

I A ( 1 )  *  IND2R (1,1) 

IA (2)  *  IND2R (2,1) 

K  -  3 

IF  (Rl  .EQ.  1)  GOTO  60 
DO  50  I  *  1 ,R1 
DO  40  J  *  1,2 
M  =  K  -  1 
DO  30  L  ■  1,H 

IF  ( IND2R ( J , I )  .EQ.  IA(L>>  GOTO  40 
30  CONTINUE 

I A (K)  *  IND2R ( J , I ) 

K  =  K  +  1 
40  CONTINUE 

50  CONTINUE 
60  K  ■  K  -  1 


n  n 


90 


100 

no 

120 


...  Remake  the  reduced  matrix  IND1  -  for  each  column 
in  IND1  (corresponding  to  a  data  set)  the 
entries  are  zero  except  for  the  entries  corres¬ 
ponding  to  the  left  array  <  =  1)  and  the  right 
array  <=  2) . 

DO  80  I  =  1,R1 
DO  70  J  =  1,K 

IND1 <J,I>  =  0 

IF  ( IND2R (1,1)  .EQ.  IA(J>)  IND1(J,I)  =  1 
IF  ( l ND2R (2,1)  .EQ.  IA<J>>  IND1<J,I>  =  2 
CONTINUE 
CONTINUE 

...  Reduce  the  arrays  CROSSA  and  MEAN  to  account 
for  the  removed  data  sets. 

DO  120  I  *  1 , R1 
DO  90  J  =  1,6 

MEAN ( I  ,  J  )  =  MEAN (DATSET < I > ,J ) 

CONTINUE 
DO  llO  J  =  1,3 
DO  100  L  =  1,3 

CROSSA (I , J,L>  «  CROSSA ( DATSET ( I ) , J , L ) 
CONTINUE 
CONTINUE 
CONTINUE 
RETURN 
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